Non-axisymmetric instability and fragmentation of general relativistic quasi-toroidal stars 
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In a recent publication, we have demonstrated that differentially rotating stars admit new chan- 
nels of black hole formation via fragmentation instabilities. Since a higher order instability of this 
kind could potentially transform a differentially rotating supermassive star into a multiple black 
hole system embedded in a massive accretion disk, we investigate the dependence of the instabil- 
ity on parameters of the equilibrium model. We find that many of the models constructed exhibit 
non-axisymmetric instabilities with corotation points, even for low values of T / 1 W|, which lead to a 
fission of the stars into one, two or three fragments, depending on the initial perturbation. 

At least in the models selected here, anm = 1 mode becomes unstable at lower values of T/\W\, 
which would seem to favor a scenario where one black hole with a massive accretion disk forms. In 
this case, we have gained evidence that low values of compactness of the initial model can lead to a 
stabilization of the resulting fragment, thus preventing black hole formation in this scenario. 



PACS numbers: 04.40.Dg, 04.25.Dm, 97.60.Lf, 04.70.-s 
I. INTRODUCTION 

The study of oscillations and stability of stars has a 
long history (see e.g. [1]). While its classical results 
tried to address the limits of possible models for main- 
sequence stars, and the question of binary star forma- 
tion from protostellar clouds, it has been extended in the 
last century to examine the properties of relativistic fluid 
equilibria, and their connection to the formation of black 
holes by axisymmetric instabilities. 

This work is an extension of our publication 12], 
where we have studied the formation of black holes 
from fragmentation in N = 3 polytropes, which was in- 
duced by the growth of a non-axisymmetric instability 
in a strongly differentially rotating star. We would like 
to focus attention here on the question of whether the 
instability we have observed is generic for strongly dif- 
ferentially rotating polytropes, and how changes in pa- 
rameters of the initial model affect its development. In 
addition, we will present some tentative results related 
to arguments previously presented by Watts et al. [3] on 
a possible connection between low-r/| and spiral- 
arm instabilities and the location of the corotation band 
in a sequence of increasing rotational energy. 

Due to the recent prospect of detecting gravitational 
radiation directly, the connection between the local dy- 
namics of collapse and gravitational wave emission is 



Here, and subsequently, T/|W| shall denote the ratio of rotational 
kinetic to gravitational binding energy in the axisymmetric initial 
model. For a definition of these quantities, see | 4]. 



currently receiving increased attention (e.g. [S', '6^). In 
this context, a non-axisymmetric instability in a star is 
expected to change the nature of the signal, and to en- 
hance the chances of detecting it [7]. We shall discuss a 
number of scenarios for gravitational collapse and black 
hole formation to illustrate this point. 

1. Stars retaining spherical symmetry: If the initial mat- 
ter distribution has spherical symmetry, no grav- 
itational waves are emitted as a consequence of 
Birkhoff's theorem. The exact solution by Op- 
penheimer and Snyder [8] already exhibits many 
features of the local dynamics, while the connec- 
tion to dynamical stability in general relativity has 
been made explicit by Chandrasekhar [9]. The 
assumption of spherical symmetry, whilst restric- 
tive, already admits simple models of phenomena 
like mass limits for compact stars, some generic 
properties of black hole formation from supermas- 
sive stars and neutron stars (e.g. ifTol Oil Mm, 
and the dynamics of apparent and event horizons 
(ibido). 

2. Stars retaining (approximate) axisymmetry: If the 
symmetry assumption is relaxed to axisymmetry, 
models of gravitational collapse admit a num- 
ber of additional features, the most important in 
this context being the emission of gravitational 
radiation. Stars in axisymmetry can be rotat- 
ing, which changes the radial modes of the non- 
rotating member of a sequence into a quasi-radial 
mode. If that is unstable, the star may collapse to 
a black hole in a manner which is similar to the 
spherically symmetric case in its bulk properties. 
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and it proceeds by (i) contraction due to a quasi- 
radial instability, (ii) formation of an event hori- 
zon centered on the axis, and (iii) ring-down to 
a Kerr black hole with a disk. We will call this 
process the canonical scenario to represent that it 
provides the expected properties of the collapse 
of slowly rotating stars. It has been studied ex- 
tensively in numerical investi gati ons |6, 13) Jl3.[T^ 
[M [H [S H SI in MM. 1261 iMfand 
seems even appropriate to describe the quality of 
collapse of most rapidly and differentially rotating 
neutron stars [25]. It should not be considered im- 
plicit here that the canonical scenario is generic for 
axisymmetric collapse: see e.g. J2^l29|| for systems 
involving toroidal black holes. 

3. Stars not retaining (approximate) axisymmetry: As al- 
ready mentioned, even in many numerical mod- 
els with three spatial dimensions, the collapse to 
a black hole proceeds in an almost axisymmetric 
manner, although the initial data is represented on 
discrete Cartesian grids. It has been found that 
even when non-axisymmetric perturbations are 
applied to the collapsing material, no large devi- 
ations from axisymmetry are seen during the col- 
lapse [26]. Judging from the perturbative theory 
of Newtonian polytropes [1], this indicates that ei- 
ther the amount of rotational over gravitational 
binding energy r/|W| is insufficient, or that the 
collapse time is too short to admit growth of initial 
deviations from the symmetric state to significant 
levels. 

The situation can be quite different when the sys- 
tem is not unstable to axisymmetric perturbations, 
or if the collapse stabilizes around a new equilib- 
rium with higher r/|W|. The classical limit of 
r/|W| w 0.27 for Maclaurin spheroids [30] indi- 
cates the onset of a dynamical instability to tran- 
sition to the X = +1 Riemann S-type sequence 
| |30l l3l|[ . To which extent this idealized behaviour 
is also realized in general relativistic compressible 
polytropes, and, more specifically, how it is con- 
nected to the formation of black holes, is the issue 
we would like to address in part here. 

If a general relativistic star encounters a non- 
axisymmetric instability, the nature of its subsequent 
evolution may be characterizable by certain properties 
of the equilibrium model, like the rotation law, T/| W|, 
compactness, and equation of state. For the limit of 
uniformly rotating, almost homogeneous models of low 
compactness, we expect, for r/|W| > 0.27, a dynami- 
cal transition to an ellipsoid by a principle of correspon- 
dence with Newtonian gravity. 

By relaxing all but the assumption of low compact- 
ness, we can make use of the rich body of knowledge 
about the stability and evolution of stars in Newtonian 
gravity. Two classical applications of stability theory 
are the oscillations of disks and the fission problem [1]. 



These matters have been investigated extensively over 
the last decades g H HI IMM HE HI HI, IS IS 
mmiilliHil^El, and a nxmiber of possible scenar- 
ios have emerged for the non-linear evolution of non- 
axisymmetric dynamical instabilities in rotating poly- 
tropes: 

1. The polytrope develops a bar-mode instability 
similar to the Maclaurin case, and possibly retains 
this shape over many rotational periods (e.g. |'4l'|). 

2. Two spiral arms and an ellipsoidal core region de- 
velop, where the latter transports angular momen- 
tum to the spiral arms by gravitational torques, is 
spun down, and collapses (e.g. [33]). This scenario 
is interesting for black hole formation, since the 
rotational support of the initial model can be par- 
tially removed by such a mechanism. If this trans- 
port is efficient enough, the core ellipsoid might 
collapse, resulting in a Kerr black hole with a disk 
of material around it. One might conjecture that, if 
the equation of state is soft enough, the disk itself 
may be subject to fragmentation, and form sev- 
eral smaller black holes which are subsequently 
accreted. 

3. One spiral arm develops, and the mode saturates 
at some amplitude js^lisll , leaving a central con- 
densation. This might also lead to central black 
hole formation. Note that the onset for this dy- 
namical instability in terms of T/| W| can be sig- 
nificantly lower than the Maclaurin limit. 

4. For polytropes with strong differential rotation, 
the initial model may be quasi-toroidal, i.e. it has at 
least one isodensity surface which is homeomor- 
phic to a torus. If models of this kind, or purely 
toroidal ones, are subject to the development of 
a non-axisymmetric instability, they may exhibit 
fragmentation |40, 43]. This is clearly the most in- 
teresting setting for the fission problem, but has 
been discussed also in the context of core collapse 
("collapse, pursuit and plunge" scenario. Fig. 24.3 
in [49], also [37, 45]). 

It is this last kind of instability we will investigate here 
in the context of general relativity, and its relation to the 
formation of black holes. 

Concerning the nature of the spiral-arm and low- 
T/|W| instabilities. Watts et al. ((1, see also I© H) 
have suggested a possible relation of these instabilities 
to their location with respect to the corotation band^. 



^ The corotation hand in a differentially rotating star is the set of fre- 
quencies associated with modes having at least one corotation point, 
i.e. a point where the local pattern speed of the instability matches 
the local angular velocity. 
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That corotation has a bearing on instabilities in dif- 
ferentially rotating disks has been suggested for some 
time; the perturbation operator is singular at corota- 
tion points, which gives rise to a continuous spectrum 
of "modes." While the initial-value problem of per- 
turbations associated with the continuous spectrum of 
stars is not well imderstood even in Newtonian gravity, 
there is evidence |50, 52, 53, 54] that a mode entering 
corotation may be subject to a shear-type instability, or 
that it merges with another mode inside the corotation 
band, which appears to admit a certain class of solu- 
tions showing similar properties as the solutions in the 
discrete spectrum [51J. While three-dimensional New- 
tonian simulations would appear, at least as a first step, 
most appropriate to gain more intuition in these matters 
155], we will collect some evidence on corotation points 
in the evolutions presented here as well. 

Since the parameter space of possible initial models 
is large, and given that three-dimensional simulations 
of this kind are still quite expensive in terms of com- 
putational resources, we restrict attention to several iso- 
lated sequences, where just one initial model parameter 
is varied to gain evidence on its systematic effects, and 
to a plane in parameter space defined by a constant cen- 
tral rest-mass density and a fixed parameter F = 4/3 
in the F-law equation of state P = (F — l)pe. We will 
find that, at least as long as we are concerned with the 
question when certain modes become dynamically un- 
stable on a sequence, the consideration of models of con- 
stant central density is not overly restrictive as far as 
the development of the instability is concerned, while 
the nature of the final remnant might be rather sensi- 
tive to it. This latter issue, namely whether a black hole 
forms or not, will not be answered in full here, since 
we will only determine whether the fragment stabilizes 
during collapse, and re-expands, or if it does not. We 
leave the location of the apparent horizon with adaptive 
mesh-refinement techniques and the subsequent evolu- 
tion with excision to future work, and concentrate here 
on the general structure of the parameter space and its 
relation to the non-axisymmetric instability. 

The choice F = 4/3 is well known to approximately 
correspond to the adiabatic coefficient of a degenerate, 
relativistic Fermi gas or a radiation pressure dominated 
gas [56], and is thus closely connected to iron cores and 
supermassive stars. We would like to point out that a 
collapsing iron core, even if its initial state is assumed 
to be determined mostly by electron pressure, is subject 
to a complex set of nuclear reactions, which involve the 
generation of neutrinos and transition to nuclear matter 
at high densities. It is because of this complexity, which 
we do not take into account here in order to reduce the 
number of free parameters, that we do not suggest the 
use of the fragmentation and black hole formation in- 
vestigated in [2] as a highly idealized model of core col- 
lapse. For supermassive stars, the situation is different, 
since an event horizon can form before thermonuclear 
reactions become important, depending on the metallic- 



ity and mass of the progenitor fsT]. Because of this, it is 
conceivable that the type of evolution in [2] can be used 
as an approximate model of supermassive black hole 
formation. Finally, we would like to mention the pos- 
sibility that gravitational wave detection may uncover 
so-far imexpected processes involving black hole forma- 
tion, and in that case it is useful to have a general under- 
standing of possible dynamical scenarios. 

II. PREVIOUS WORK 

The background for this study comes from three ar- 
eas: the study of (i) fragmentation in Newtonian poly- 
tropes, (ii) non-axisymmetric instabilities in general rel- 
ativistic polytropes, and (iii) black hole formation by 
gravitational collapse. The first area is represented by 
a large number of publications (some of them cited in 
the introduction), but we would like to mention specifi- 
cally the work by Centrella, New et al. (43l|58j], since the 
kind of initial model and subsequent evolution studied 
in these publications are similar to the ones presented 
by us in |j2|] and here, apart from the fact that Newto- 
nian gravity and a softer equation of state (F = 1.3) 
was used. New and Shapiro [59] investigated equi- 
librium sequences of differentially rotating Newtonian 
polytropes with F = 4/3 to present an evolutionary 
scenario where supermassive stars develop a bar-mode 
instability instead of collapsing axisymmetrically. This 
kind of scenario (see also [2/,60J) is also important when 
connecting the general relativistic fragmentation pre- 
sented here to the evolution of supermassive stars. 

Non-axisymmetric dynamical instabilities in general 
relativistic, self-gravitating fluid stars have been stud- 
ied by several authors [2al6ll|6l,[6l[6l|65*|. Some evi- 
dence of fragmentation has been found in [25] in a ring 
resulting from a "supra-Kerr" collapse with } /M^ > 1 
(here / denotes the total angular momentum, and M 
the ADM mass), but no black hole was identified, al- 
though the pressure in the initial data was artificially re- 
duced by a large factor in order to induce collapse. Fi- 
nally, black hole formation by gravitational collapse has 
been studied extensively, (see references in the introduc- 
tion], and in recent years also in three spatial dimensions 
lHllIIIllllllIEIllllEI^- The collapse of differen- 
tially rotating supermassive stars in the approximation 
of spatial conformal flatness has been investigated by 
Saijo (H. 

In addition, the work on low-T/|W| instabilities by 
Watts et al. [3, 50, 51] has already been described in the 
introduction, and recent numerical studies of related in- 
terest can be found in L5^[6^[69ll . 

III. NUMERICAL METHODS 

All model calculations presented here are performed 
in full general relativity. We use Cartesian meshes 
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with the only symmetry assumption being reflection 
invariance with respect to the equatorial plane of our 
models. We employ the Cactus computational frame- 
work fTO", 'tT] in our simulations. We use finite- 
differencing methods of second order for the spacetime 
variables, and finite-volume, high-resolution shock- 
capturing techniques for the hydrodynamical variables. 
Hydrodynamics and metric evolution are coupled by 
means of the method of lines, using the time-explicit itera- 
tive Crank-Nicholson integrator with three intermediate 
steps IzIIzSH. 

In the following we give a brief overview on the 3- 
metric evolution, hydrodynamics and mesh refinement 
methods used. We set c = G = K = 1 to fix our system 
of units, where c denotes the speed of light, G the grav- 
itational constant and K the polytropic constant in the 
relation P = Kp^ between pressure P and rest-mass den- 
sity p used to compute the initial models. Latin indices 
run from 1 to 3, Greek from to 3. We use a spacelike 
signature (—,+,+,+). Unless explicitly mentioned oth- 
erwise, we assume Einstein's summation convention. 



A. Metric evolution 

We evolve the 3-metric with the 3+1 Cauchy free- 
evolution code developed at the Albert Einstein Insti- 
tute f/^, '75']. The code provides at each time step a so- 
lution to the Einstein equations in the ADM 3+1 formu- 
lation [76], while internally evolving the spacetime in 
the NOK-BSSN conformal-traceless recast of the ADM 
equations [T^ [TSL 79] which has been proven to lead 
to much more stable numerical evolutions of Einstein's 
equation than the original ADM system [72]. The de- 
tails of our NOK-BSSN implementation can be found in 
ll72l l801. Here we only briefly summarize the ADM sys- 
tem through which we couple spacetime and hydrody- 
namics II271I . 

In ADM, the four-dimensional spacetime is foliated 
into a set of three-dimensional and non-intersecting 
spatial hypersurfaces. Individual surfaces are related 
through the lapse function a, which describes the rate 
of advance of time along a timelike unit normal vec- 
tor n*, and through the shift vector /5' which indicates 
how the coordinates change from one slice to the next. 
The gauge-freedom in ADM allows a free choice of both 
lapse and shift, but care must be taken for the choice of 
gauge may have an effect on numerical stability [8Q1 . 

The ADM line element reads 



^i^')df + 2^idx'dt + jjjdx'dx! (1) 



and the first-order form of the evolution equations for 
the spatial 3-metric Jij and the extrinsic curvature Kjj 
read 



dtjij = -2t^K,j + \/,^j + \/jlii, 



Rij + KK,j 



(3) 



where V, denotes the covariant derivative with respect 
to the coordinate direction 3,- and the 3-metric 7/y. Rjj is 
the coordinate representation of the 3-Ricci tensor and 
K = Y^^ij is the trace of the extrinsic curvature. Padm is 
the energy density as measured by an Eulerian observer 



I8III . Sij is the projection of the stress-energy tensor on 
the spacelike hypersurfaces and S = Y^S 



•J- 



Gauge choices 

We evolve the lapse function with the hyperbolic K- 
driver condition I80ll82i1 of the form 



dtoc 



-a^f{ci){K-Ko), 



(4) 



where /(a) is an arbitrary positive function of a and 
Kq = K{t = 0). We choose /(a) =11 ol which is referred 
to as l+/o^ slicing and has excellent singularity-avoiding 
properties in the sense that the lapse tends to zero near a 
physical singularity, freezing the evolution in that region. 

For the model simulations presented in this paper we 
find that a dynamical evolution of the spatial gauge /5' is 
not necessary, and we keep it fixed to its initial direction 
and magnitude. 



B. General-relativistic hydrodynamics 

The equations of general-relativistic hydrodynamics 
are derived from the conservation equations of the 
stress-energy tensor T"^ and the matter current density 
f: 



v„r'' = o, v„f = o. 



(5) 



where /" = pu", p is the rest-mass density and u" the 
4-velocity of the fluid. We use the perfect-fluid stress 
energy tensor 



rab 



phu'ii'' + Pg' 



ab 



(6) 



(2) 



with P being the fluid pressure, h = 1 + € + P / p the 
relativistic specific enthalpy and e the specific internal 
energy of the fluid. 

For evolving the hydrod3mamical fields we employ 
the Whisky code [27, 83] which implements the general 
relativistic hydrodynamics equations in the hyperbolic 
first-order flux-conservative form proposed and tested 
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in [84, "S^. This code requires us to add an artificial 
atmosphere to the computational domain in regions of 
very low density. We typically choose an atmospheric 
density of 10~^ of the maximal density of the initial 
model. The evolved state vector = (D, S„ t)^ is de- 
fined in terms of the primitive hydrodynamical variables 
p, P, and v', the 3-velocity, measured by an Eulerian ob- 
server: 



u = 


' D ' 










T 





^phW^Vj 
^{phW^ - P - Wp) 



DiU" 



where 7 = det^ij and W 

Lorentz factor. 

The set of equations then reads 



(7) 



1/1/1 — jijV^v! is the 



dtU + 8,F' = S, 
with the three flux vectors given by 

a{{v'-l[i')Sj + ^P3p 



(8) 



(9) 



The source vector S, which contains the curvature- 
related force and work terms, but no derivatives of the 
primitive variables, is given by 







(10) 



where V^fw are the standard 4-Christoffel symbols. 
We choose the ideal fluid T-law equation of state, 

Pip,e) = {T-l)pe (11) 

to close the system of hydrodynamic equations. 

C. Mesh refinement 

In order to ensure adequate spatial resolution whilst 
keeping the computational resource requirements of our 
three-dimensional simulations to a minimum, we use 
Berger-Oliger style l86ll mesh-refinement with subcy- 
cling in time as implemented by the open-source Carpet 
iSZ, 88] driver for the Cactus code. Carpet provides 
fixed, progressive |6] and adaptive mesh refinement. 
In this study we use a predefined refinement hierarchy 
with five levels of refinement, arranged in a box-in-box 
manner centered on the origin. The resolution factor be- 
tween levels is two. We point out that adaptive (or at 
least progressive) mesh refinement will be necessary to 
track black-hole formation in detail, and has been per- 
formed on one model in [2]. 



D. Mode extraction 

To evaluate and quantify the stability or instability of 
a given model to non-axisymmetric perturbations, we 
extract azimuthal modes e™*' by means of a Fourier 
analysis of the rest-mass density on a ring of fixed co- 
ordinate radius in the equatorial plane''. Following 
Tohline et al. t32ll , we compute complex weighted av- 
erages 



C„ 



2n 



p{cD, cp,z = 0)e''"1'd(p 



2tt Jo 

and define normalized real mode amplitudes 

Cm 



An 



Co 



(12) 



(13) 



Here cD — + = const, and is chosen to corre- 
spond to the initial equatorial radius of maximum den- 
sity in our quasi-toroidal models, if not mentioned oth- 
erwise. The index m corresponds to the number of az- 
imuthal density nodes and is used to characterize the 
modes. 



IV. INITIAL DATA 

A. Quasi-toroidal polytropes 

We focus on general relativistic, differentially rotating 
polytropes which are quasi-toroidal: Such a polytrope has 
at least one isodensity surface which is homeomorphic 
to a torus. To construct equilibrium polytropes of this 
kind, an extended version of the Stergioulas-Friedman 
(RNS) code is used [4, 89, 90], which uses the numerical 
methods developed in Komatsu, Hachisu and Eriguchi 
ll9lll93.l93ll . The code assumes a certain gauge in a sta- 
tionary, axisymmetric spacetime, such that we can write 
the line element in terms of potentials v, ip, to and }i, and 
the Killing fields dt and as [4] 



ds^ 



e^^dt^ + e^V{d(p - oodtf + e^f{dr^ + r^dO^) 



(14) 

To compute an equilibrium polytrope, the central rest- 
mass density pc, the coordinate axis ratio fp / and a 
barotropic relation P{p) need to be specified^. In addi- 
tion, the equations of structure [4] contain an additional 
freely specifiable function F{Cl). We will use the com- 
mon choice 



F(n) 



(15) 



^ These quantities are not gauge-invariant, but they provide a useful 
way of characterizing the representation of the instability within our 
choice of coordinates. 

* Purely toroidal models have Vp / = 0. 
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where Cic denotes the angular velocity at the star's cen- 
ter. This rotation law reduces to uniform rotation in the 
limit A ^ 00, and to the constant specific angular mo- 
mentum case in the limit A ^ 0. We will, however, 
often use the normalized parameter A = A/Vf,, where 
is the coordinate radius of the intersection of the stellar 
surface with the equatorial plane = tz/2. For construc- 
tion of a polytrope, the equation of state is constrained 
to the polytropic relation 



(16) 



with the polytropic constant K and the coefficient T, 
which can also be expressed by the polytropic index 
N = (F — Without loss of generality we set K = 1 
in all cases. The equation of state of the initial model is 
completed by the energy relation e = P/[(F — l)jD]. 

The resulting set of equations is solved iteratively |4], 
where the initial trial fields are a suitable solution of the 
TOV system. To converge to the desired model, it may 
be necessary to select a number of intermediate attrac- 
tors as trial fields. Some models are thus constructed by 
first obtaining a specific quasi-toroidal model, and then 
moving in parameter space in the quasi-toroidal branch 
to the target model. In this work, a hook model with pa- 
rameters = 0.3 and {rp/ri,)hook = 0-15 is gener- 
ated, which then is used as initial guess to construct the 
target model. 

If we include the polytropic coefficient F, we 
have to consider a four-dimensional parameter space 
{T ,pc, A,rp/re). We will not study the whole parame- 
ter space here: Rather, we first use a reference model 
and explore sequences in pc, F and Yp / containing this 
model, and will subsequently concentrate on the impor- 
tant case F = 4/3, since it approximately represents a 
radiation-pressure dominated star. 

Most polytropes have been constructed with a merid- 
ional grid resolution of n,- = 601 radial zones and 
"cos 9 = 301 angular zones, a maximal harmonic index 
f-max = 10 for the angular expansion of the Green func- 
tion and a solution accuracy of lO^"^. Selected models 
have been tested for convergence with resolutions up to 
rir = 2401, «cose = 1201, and (.,nax = 20. 

To investigate the stability of the polytropes con- 
structed with the RNS code, two kinds of perturbations 
are applied: the pressure is reduced by 0.1%, and a cylin- 
drical density perturbation of the form 



p{x)^ p{x) 



re 



\n,Bf{cD)sm{mcp) (17) 



m=l 



is added to the equilibrium polytrope. Here, m e 
{1,2,3,4}, A m is either or 1, is the cylindrical radius 
and f{cD) is a radial trial function. Experiments have 
been made with f{cD) = cD and f{cD) = CD"', but the 
exact choice was found not to affect the results signifi- 
cantly. This is true quite generally, since we only require 
the trial fimction to have some reasonable overlap with 
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Figure 1: Time evolution of the norm of the Hamiltonian 
constraint for different resolutions (the numbers refer to grid 
point on a single patch of our mesh refinement grid). The time 
is normalized to the dynamical timescale = Re\/Re/M. 



a set of quasi-normal modes. It is beyond our scope to 
investigate the full spectrum of quasi-normal modes of 
general relativistic polytropes; therefore, we determine 
stability only with respect to specific trial functions. The 
choice made is not completely arbitrary, however: a 
quasi-toroidal polytrope has an off-center toroidal re- 
gion of maximal density, and it is this region which 
will dominate the dynamics if a fragmentation instabil- 
ity sets in. A linear perturbation without nodes in this 
region can be expected to be compatible (have non-zero 
scalar product) with most low-frequency quasi-normal 
modes. The function /(co) = CD'" has the additional 
property of smoothness at the center, but, as already 
noted, numerical experiments have shown the differ- 
ence to be negligible in practice. An additional note on 
the use of language: If we find that a perturbation with 
A; = Sjj, i E {1,2,3,4} leads to an instability with the 
associated number of node lines in the equatorial plane, 
we will denote this instability with the term m = j mode 
(and the corresponding perturbation m = j perturbation). 
This is a simplification insofar as each m is expected to 
represent a (discrete) infinite spectrum of modes j|93l , 
from which we will observe only the fastest-growing 
unstable member. While we will attempt to discuss the 
nature of the global evolution to some extent, we will, 
for the reasons stated above, concentrate on the high- 
density torus, and mostly neglect the dynamics of its 
halo. 

The perturbations applied are both constraint- 
violating. This is no significant issue, since the discrete 
evolution of the NOK-BSSN system will introduce con- 
straint violations even if some minimization technique 
has been applied to the initial data. Fig.[l]shows the evo- 
lution of the norm of the Hamiltonian constraint for 
different grid resolutions, for a model perturbed with 
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Figure 2: Decadic logarithm of the der\sity ir\ a meridional 
plane of the model constructed with the parameters in TableU 
The model is of quasi-toroidal nature. Note that this contour 
plot is included to demonstrate the qualitative properties of 
the initial model. Therefore, no detailed scale is provided. 



A, = Sii and B = 10 ^. (The time in this figure, and in 
the subsequent ones, is normalized to t^y,, = RcvRT/M, 
where Rf, is the circumferential equatorial radius, and 
M is the ADM mass.) Note that, for a typical pertur- 
bation amplitude of Sp/ p w 10"'', the Hamiltonian con- 
straint will be violated by 3H w 1671,0 ■ 10"^ w 2.5 ■ 10"^, 
which is sufficiently smaller than the violations during 
evolution in Fig. [T] Also, tests have been performed 
where the perturbation amplitude B (cf . eqn. [T7)l is re- 
duced by a factor of 10, and found that this does not 
affect the growth rate of the perturbation, as expected 
from small perturbations. To conveniently compare dif- 
ferent resolutions, the amplitude is kept constant; how- 
ever, it is possible to reduce the perturbation ampli- 
tudes for resolutions significantly higher than the ones 
used here (specifically, in the regime where SH would 
become comparable to the constraint violations during 
evolution) to obtain a system where convergence, now 
to the equilibrium system, does only depend on the 
well-posedness of the continuum initial-boundary value 
problem and the stability of the discrete system^. 



B. The reference polytrope and associated sequences 

We start with a polytrope with the same central rest- 
mass density (pc = 3.38 ■ 10~^) as Saijo's series of dif- 
ferentially rotating supermassive star models [26!]. To 
obtain experience with the influence of certain parame- 
ters on the stability properties of the relativistic quasi- 
toroidal polytropes, some sequences containing this ref- 
erence model have been constructed. 



^ In addition, since our technique of solving the equations of hydro- 
dynamics require us to add an artificial atmosphere in the vacuum 
region, one would also need to reduce its density with resolution. 



Polytropic index T 


r 


4/3 


Central rest-mass density 


Pc 


3.38 ■ lO-*^ 


Degree of differential rotation 


A 


1/3 


Coordinate axis ratio 


rp/re 


0.24 


Density ratio 


Pmax/ Pc 


16.71 


ADM mass 


M 


7.003 


Rest mass 


Mo 


7.052 


Equatorial inverse compactness 


Re/M 


11.71 


Angular momentum 


J 


52.20 


Normalized angular momentum 


J/M^ 


1.064 


Kinetic over binding energy 


T/\W\ 


0.227 


(See caption) 


Of. /O^ 


0.467 



Table I: Parameters and integral quantities of the reference 
quasi-toroidal polytrope |2]. The quantities T, pc, A and Vp/rp 
are parameters. The quantity Cle is the angular velocity on 
the equator, while is the associated Keplerian velocity of 
the same model. Therefore, the mass-shedding sequence is lo- 
cated atne/Qjc = 1. 



1. The reference model 

The reference model is identical to the model used in 
[2]; its parameters and integral quantities are shown in 
Table U Fig.|2]is a graph of the density in the meridional 
plane for this model. This model has been found to be 
unstable to perturbations with m = 1 and m = 2 (see be- 
low, and [2]), which lead to fragmentation. In the case of 
m = 1, black hole formation has been demonstrated by 
locating an apparent horizon centered on the fragment 

m. 



2. The sequence of axis ratios 

A number of sequences containing this model have 
been constructed to study the stability properties when 
varying typical parameters. The first of these sequences 
is parameterized by the coordinate axis ratio rp/r^. Its 
members are denoted by R<rp / re>, and their properties 
can be found in Table [III It is apparent that, with in- 
creasing axis ratios, the quantity T/ 1 W | and the ratio of 
maximal to central rest-mass density Pmax/pc both de- 
crease monotonically. None of the sequence members is 
close to the mass-shedding limit. Below rp/r^ = 0.20, 
no models could be constructed due to failure of con- 
vergence. We will also not consider models with larger 
Ypl Yf.. The reason for this is indicated in Fig.|3l Beyond 
an axis ratio of rp/r^ w 0.3423, quasi-toroidal models 
could not be constructed as the numerical code fails to 
converge. 
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Model 


P 


c 




rp/re 


pmax/ Pc 


M 


Mo 


Re/M 


/ 


//M2 


T/\W\ 




R0.20 


3.38- 


10- 


-6 


0.20 


38.12 


6.181 


6.200 


9.660 


38.59 


1.010 


0.235 


0.487 


R0.22 


3.38- 


10- 


-6 


0.22 


25.69 


6.662 


6.710 


10.41 


45.46 


1.024 


0.228 


0.475 


ROM 


3.38- 


10- 


-6 


0.24 


16.76 


6.989 


7.037 


11.71 


52.00 


1.065 


0.227 


0.467 


R0.26 


3.38- 


10- 


-6 


0.26 


11.07 


7.334 


7.391 


13.10 


58.99 


1.097 


0.223 


0.460 


R0.28 


3.38- 


10- 


-6 


0.28 


7.312 


7.585 


7.646 


14.83 


65.19 


1.133 


0.219 


0.455 


R0.30 


3.38- 


10- 


-6 


0.30 


4.733 


7.764 


7.825 


17.13 


70.82 


1.175 


0.213 


0.452 


R0.32 


3.38- 


10- 


-6 


0.32 


2.934 


7.847 


7.905 


20.46 


75.48 


1.226 


0.207 


0.452 


R0.34 


3.38- 


10- 


-6 


0.34 


1.539 


7.755 


7.803 


27.42 


78.72 


1.309 


0.196 


0.463 



Table II: Parameters and mtegral quantities of the R sequence of axis ratios, which contains the reference model for rp/r^ = 0.24. 
Each member of the sequence is denoted by the term R<rp/re>. All models have F = 4/3 and A = 1/3. 



Model 


r 


Pc 


pmax/ Pc 


M 


Mo 


Re/M 


/ 


J/M^ 


T/\W\ 


a/n^c 


G1.333 


4/3 


3.38-10-^ 


16.76 


6.989 


7.037 


11.71 


52.00 


1.065 


0.227 


0.467 


G1.4 


1.4 


1.32 ■ 10-5 


16.44 


2.805 


2.856 


11.69 


8.337 


1.059 


0.211 


0.424 


G1.45 


1.45 


3.7-10-5 


13.83 


1.624 


1.662 


11.70 


2.767 


1.049 


0.202 


0.404 


G1.5 


1.5 


9.2 ■ 10-5 


11.53 


1.038 


1.066 


11.65 


1.115 


1.035 


0.194 


0.390 


G1.6 


1.6 


3.75 ■ 10-* 


8.349 


0.5197 


0.5363 


11.80 


0.2758 


1.021 


0.183 


0.370 


G2.0 


2.0 


8.2 • 10-3 


3.880 


0.1270 


0.1323 


11.69 


0.0155 


0.9633 


0.159 


0.331 



Table III: Parameters and integral quantities of the G sequence of polytropic coefficient F, which contains the reference model for 
F = 4/3. Each member of the sequence is denoted by the term G<Y>. The central density is adjusted to yield approximately the 
same inverse equatorial compactness Re/ M as in the reference model. The models G1.7 to G1.9 have not been constructed, since 
the models G1.6 and G2.0 were found to be stable. All models have A = 1/3 and rp/r^ = 0.24. 



3. The sequence of stiffnesses 

This sequence is a variation of the parameter F in the 
polytropic relation P = Kp^ , which also determines 
the stiffness of the ideal fluid equation of state P = 
(r — l)pe used for evolution. To obtain a sequence of 
comparable compactness, we adjust the central density 
Pc to yield approximately the same Re/M. The parame- 
ters and integral quantities are shown in Tablelllll Along 
the sequence of increasing T, the value T/| W| decreases 
from 0.227 to 0.159. Therefore, it would also be inter- 
esting to consider a sequence of models with varying T, 
but constant T / 1 W | , by adjusting the axis ratio / ac- 
cordingly. Unfortunately, such a sequence could not be 
obtained, since the initial data solver did not converge 
to models with the required (low) axis ratios. 

4. The sequence of compactnesses 

The next sequence is a variation of the central rest- 
mass density pc while leaving all other parameters fixed. 
The resulting sequence, as is apparent from Table HVl is 
also a sequence of inverse equatorial compactnesses Re/M. 
The members have been selected to represent models 
which are half to one eighth as compact as the refer- 
ence polytrope. The sequence shows that T/| W| is only 
slightly affected by the choice of pc, but Re/M and // 



change significantly. 



C. Quasi-toroidal and spheroidal models of constant 
central rest-mass density 

In addition to sequences containing the reference 
model, we explore a more extended part of the param- 
eter space of models. We use T = 4/3 and pc = 10^^ 
to define a surface in the parameter space spanned by 
the axis ratio rp / re and the degree of differential rota- 
tion A. While an ideal fluid with T = 4/3 is an approxi- 
mation for the material properties of radiation-pressure 
dominated stars, the restriction to pc = 10^'' is arbi- 
trary. However, as discussed in Section |Vl the non- 
axisymmetric stability of the quasi-toroidal models is 
probably less sensitive to pc than to rp/re or A. The re- 
striction to a plane is necessary since three-dimensional 
simulations of quasi-toroidal relativistic stars are still ex- 
pensive; however, selected models will also be studied 
with different pc in Section FVl 

The general properties of the polytropes obtained 
with the RNS code is shown for a central rest-mass den- 
sity pc = 10^^ in Fig. The top left plot shows the 



* We have also generated the same set of plots for the plane pc = 
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Model 


Pc 


rp/re 


pniaxf Pc 


M 


Mo 


Re/M 


/ 


//M2 


T/|W| 




CI 


3.38 ■ 


0.24 


16.76 


6.989 


7.037 


11.71 


52.00 


1.065 


0.227 


0.467 


C2 


1 ■ 10-^ 


0.24 


31.06 


10.65 


10.74 


22.72 


167.0 


1.474 


0.225 


0.434 


C4 


7.5 ■ 10-"^ 


0.24 


37.01 


12.54 


12.60 


45.65 


326.0 


2.073 


0.225 


0.423 


C8 


8-10-1° 


0.24 


39.64 


13.47 


13.51 


89.91 


525.6 


2.897 


0.225 


0.419 



Table IV: Parameters and integral quantities of the C sequence of compactnesses, which contains the reference model for pc = 
3.38 ■ 10-^. Each member of the sequence is denoted by the term C<a>, where a denotes the approximate ratio of Re/M to 
(Re/M)ref of the reference model CI. The sequence is obtained by varying the central rest-mass density. All models have 
r = 4/3 and A = 1/3. 



Model 


Pc 


A 


rp/re 


Pmaxf Pc 


M 


Mo 


Re/M 


/ 


//M2 


T/|W| 


ae/ClK 


A0.1R0.15 


10- 


7 


0.1 


0.15 


246.8 


4.896 


4.893 


20.82 


18.62 


0.777 


0.124 


0.195 


A0.1R0.50 


10- 


7 


0.1 


0.5 


1.881 


5.387 


5.390 


106.2 


31.69 


1.092 


0.0706 


0.126 


A0.3R0.15 


10- 


7 


0.3 


0.15 


356.9 


7.964 


8.034 


12.04 


70.66 


1.114 


0.228 


0.434 


A0.3R0.50 


10- 


7 


0.3 


0.5 


1.00005 


6.291 


6.300 


151.4 


77.10 


1.948 


0.108 


0.360 


A0.6R0.15 


10- 


7 


0.6 


0.15 


541.2 


21.29 


21.39 


61.62 


1360 


3.000 


0.276 


0.650 



Table V: Parameters and integral quantities of selected quasi-toroidal models in the parameter space plane defined by pc = 10 
and r = 4/3. The models are labelled by A<A>R<rp/re>. All models have F = 4/3. 



function Qe/Qj^, where is the equatorial stellar an- 
gular velocity, and Qj^ is the corresponding Keplerian 
angular velocity. The jump indicated in the equilibrium 
model surface has already been f oimd in the R sequence, 
see Fig.|3l 

The topological nature of the polytropes is shown 
in the top right panel of Fig. HI which plots the ratio 
of maximal to central rest-mass density Pmax/ Pc- This 
value measures the degree of toroidal deformation of 
the model, with the limiting Ca.SGS Pmax /pc = I (purely 
spheroidal polytrope) and pmax/ pc = ^ (purely toroidal 
polytrope). Since we are interested in the properties of 
quasi-toroidal models, we will concentrate our study on 
the part of this plot covered by contour lines. 

Judging from the study of Newtonian polytropes, one 
would expect that the function T / 1 W| is related to non- 
axisymmetric stability. For the sequence of Maclaurin 
spheroids, the dynamically unstable subset can be de- 
scribed by the simple inequality T/|W| > (T/|W|)^y„ 
| |30|] . suggesting to use r/|W| to parameterize the se- 
quence. While the situation is clearly more complicated 
with relativistic, differentially rotating polytropic mod- 
els, the middle left plot in Fig. |4] suggests that the quasi- 
toroidal models with small axis ratio rp/r^ = 0.15 are 
more likely to be unstable to non-axisymmetric pertur- 
bations. We will study this in SectionlVl 

The bottom right panel in Fig.|4]is an illustration of the 
initial model parameter space. The polytropes which 



10 and find the dimensionless quantities Cic/Cifc, pmax/ Pc, and 
T / 1 W| to be quite similar. 



have been evolved numerically are marked by circles. 
The equilibrium parameters and associated quantities of 
a selected set of these polytropes are listed in Table |Vl 



D. A sequence of central rest-mass densities containing 
the model A0.2R0A0 

The L sequence (see Table IVB is a variation of the 
C sequence in Section FlV B 41 It starts from the model 
A0.2R0.40 instead of the reference model. In contrast to 
the C sequence, we do not keep the axis ratio fp / fixed 
while varying the central rest-mass density, but rather 
the quantity r/|W|. This sequence is constructed to 
study the influence of the compactness on a model near 
the boundary to the region denoted by "I" in Fig.|25l(see 
also Section |VH). 



V. RESULTS 

The models constructed in Section |IV] have been 
evolved numerically to study their stability properties. 
We will start with discussing the reference model, and 
show that it is unstable to a non-axisymmetric pertur- 
bation which leads to black hole formation. Then, the 
sequences of axis ratios, compactness and stiffness from 
Section|IV]are studied. Finally, the parameter plane con- 
strained by F = 4/3 and pc = IQ-'^ is sampled, and the 
coordinate location of the corotation point on a sequence 
in this plane is investigated. 
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Model 


r 


Pc 


A 


rp/Ve 


pmax/ pc 


M 


Mo 


Re/M 


/ 


//M2 


T/\W\ 


a/fiK 


LI 


4/3 


10- 


4 


0.2 


0.354 


1.874 


3.621 


3.578 


14.85 


9.952 


0.760 


0.144 


0.301 


L2 


4/3 


10- 


5 


0.2 


0.378 


2.434 


5.173 


5.177 


20.67 


24.06 


0.900 


0.144 


0.280 


L3 


4/3 


10- 


6 


0.2 


0.392 


2.644 


6.505 


6.524 


35.69 


49.51 


1.170 


0.144 


0.269 


L4 


4/3 


10- 


7 


0.2 


0.4 


2.689 


7.348 


7.362 


69.41 


87.33 


1.617 


0.144 


0.264 



Table VI: Parameters and integral quantities of the L sequence. This sequence is constructed by starting from the model A0.2R0.40 
(cf. Section lTVCl , which is identical to the model L4, and keeping T/\W\ fixed while reducing the central density. 




0.20 0.22 0.24 0.26 0.28 0..W 0.32 0.34 0.36 
0.25 




0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 



r /t 

p t 



Figure 3: pmax/pc (top) and T/\W\ (bottom) for the R se- 
quence. Beyond an axis ratio of rp/r^ ^ 0.3423, the initial data 
solver converges to spheroidal models. 



A. Evolution of the reference model 

The main results from evolving the reference poly- 
trope defined in Section IIV B 21 have already been dis- 
cussed in [2]. When subject to a perturbation of the form 
given in eqn.[l71 the torus transforms into one {m = 1) or 
two {m = 2) fragments. In the case of the m — 1 pertur- 
bation, it has been shown that the fragment is partially 



covered by an apparent horizon, indicating black hole 
formation. 

In this section, we will take a more detailed look at 
this model, and discuss some technical issues relevant 
to the parameter space study below. 



1 . Development of the instability 

Fig. |5] shows the development of the non- 
axisymmetric instability in the equatorial plane when 
using a perturbation of the form given by eqn. [T7| and 
A„, = 1 for m = 1, . . . , 4. The density perturbation 
is not apparent in the initial model, but, after a few 
dynamical timescales, an instability has developed 
which entirely destroys the structure of the star. (Note 
that the artifacts outside the stellar surface are caused 
by the introduction of an artifical atmosphere.) In this 
case, one collapsing off-center fragment forms in the 
system. Judging from Fig. |5] there is a "collapse of 
the lapse," which is a well-known effect when using 
singularity-avoiding slicings, and which indicates the 
development of a black hole. 

To investigate the instability more closely, the Fourier 
mode extraction discussed in Section |lll] has been ap- 
plied to the coordinate radius of highest density in the 
initial model, which is at r = 0.25Re. We concentrate on 
this radius for reasons already discussed: different ex- 
traction radii will be considered in Section FV A 21 The 
top left panel of Fig. [6] displays the evolution of the am- 
plitude of the first four Fourier modes A„i, m = 1, . . . , 4, 
in this evolution. Although all four modes have been 
injected with the same amplitude (~ 10^*), the m = 4 
mode displays a significant initial growth of about an or- 
der of magnitude, and then oscillates around this level 
until 1 1 tijy„ ~ 6, where non-linear effects become im- 
portant. The high level of m = 4 noise is very likely an 
artefact of the Cartesian grids used for the simulation. 
Support for this argument can be obtained by consider- 
ing that the equatorial section of the grid has a discrete 
C4 symmetry, and by comparing Fig. 3 and 4 in [43], 
where results from the development of a similar non- 
axisymmetric instability, though in Newtonian gravity, 
were achieved on cylindrical and Cartesian grids. We 
note that the discrete model appears to be stable against 
this perturbation, and also against the m = 3 perturba- 
tion. The remaining modes are unstable, and the growth 
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Figure 4: Polytropic models constructed with the RNS code. The plots show the parameter plane spanned by the axis ratio 
rp/ve and the differential rotation parameter A, constrained by pc = 10~^ and T = 4/3, and resolved by construction of 6400 
(80 X 80) models. Top left: contours of the equatorial stellar over Keplerian angular velocity Cie/ (the mass shedding limit is 
the isosurface Oe/O^ = 1). Top right: decadic logarithm of the ratio of maximal over central density logio{pmax/pc)- Middle left: 
rotation parameter T/| W|. Middle right: decadic logarithm of the inverse equatorial compactness Re/M. Bottom left: normalized 
angular momentum // M^. Bottom right: Selection of evolved initial models. The thick continuous line marks the jump between 
quasi-toroidal and spheroidal models apparent in the top left to bottom left plots, the thick dashed line is the mass-shedding limit, 
and the thin dashed line indicates an approximate division between spheroidal and quasi-toroidal models. The models selected 
for numerical evolution are marked by circles. 
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Figure 5: Development of the fragmentation instability in the reference polytrope. The left sequence of plots show the decadic 
logarithm of rest-mass density in the equatorial plane, the right sequence of plots shows the lapse function. The snapshots 
correspond to times 1 1 t^y„ = (top), 6.28 (middle) and 7.48 (bottom). The initial model, perturbed by eqn.[l7]with A,„ = 1 for 
m = 1, . . . , 4, develops a spiral arm instability and a collapsing fragment. The domain outside the star initially has the density of 
the artificial atmosphere {patm ~ 5 ■ 10^^"), and the artifacts outside the star are caused by interactions of the atmosphere with 
the stellar surface and the outer boundaries of the computational domain. Here, and in all evolution sequences which follow, the 
extent of the spatial coordinate domain plotted is the same in all snapshots. Also, note that the color map is adapted to the range 
of function values in each plot. 
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Figure 6: Mode amplitudes versus time, extracted at r = 0.25Rp, the radius of highest rest-mass density in the initial model, 
for different initial perturbations. The amplitude A„, is the m-th harmonic Fourier projection of the density, normalized to the 
average value. The perturbations corresponding to each plot are (cf. also eqn.[l7): Am = 1 for m = 1, ... ,4 (top left). Am = S,„i 
(top right). Am = S„i2 (middle left). Am = S,„3 (middle right). Am = S,„i (bottom left), and Am = for m = 1, . . . , 4. For details see 
text. 




Figure 7: Similar to Fig. [S] but now for a perturbation A,„ = 
(5„,i . Shown is the decadic logarithm of the density in the equa- 
torial plane. The snapshots correspond to times t/ t^y„ = 6.28 
(top), 7.11 (middle) and 7.48 (bottom). While the m = 2 mode 
is now suppressed, the qualitative evolution is similar to that 
displayed in Fig.|5l 
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Figure 8: Similar to Fig.|5] but for a perturbation Am = S,„2. 
Shown is the decadic logarithm of the density in the equatorial 
plane. The snapshots correspond to times t/t^y„ = 6.28 (top), 
7.66 (middle) and 8.85 (bottom). Two fragments develop and 
encounter a runaway instability while orbiting each other. 
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Figure 9: Comparison of mode amplitude versus time for dif- 
ferent initial perturbations. The upper panel shows the mode 
amplitude Aj, and the lower one shows Ai- Note that, in the 
case A,„ = &„^\, the m = 2 mode is dominated by non-linear 
effects from the fragmentation. 



Figure 10: Local growth times Tm{t) := (dlnA/dt)^^ (top) 
and frequencies d(p,n I dt (bottom) for modes m = 1 and m = 1, 
both taken from simulations with the corresponding initial 
perturbation. The plots denoted by "time averaged" are run- 
ning averages over t^y,,. 



times of both modes are similar. In this specific case the 
m = 1 structure dominates the late-time evolution, and 
leads to the spiral arm structure and fragmentation vis- 
ible in Fig. |5l 

The structure of the numerical noise depends on 
the grid geometry, resolution, finite difference opera- 
tors and discrete methods for treating hydrodynamics, 
the outer boundary conditions and the artificial atmo- 
sphere. Therefore, it is important to know to which 
degree the four Fourier modes are coupled during the 
evolution. Since the initial perturbation is considered to 
be "small," which holds when compared to the m = 4 
numerical noise level discussed above, we expect that 
coupling becomes important as soon as the amplitudes 
A,n get close to imity. To determine this in the reference 
model, a number of simulations have been performed 
with perturbations of the form given by eqn. [171 but 



with A,,, = S„y for different ; e {1, ... ,4} to select in- 
dividual modes, and one simulation where no pertur- 
bation is applied. The mode amplitudes in these simu- 
lations are shown in Fig. |6l and the m = 1 and m = 2 
modes are compared to the perturbation with A„j — 1 
in Fig.|9l As long as the mode development is not domi- 
nated by another mode which at a higher amplitude, e.g. 
as in the case of the m = 2 mode in the top right panel 
in Fig. [6j the growth times are comparable for different 
perturbations. 

The high-amplitude, strongly non-linear develop- 
ment at late times is sensitive to the perturbation func- 
tion, as is visible from the evolutions shown in Fig. [7| 
and |8l In the case of the m = 1 perturbation, a single 
fragment develops and collapses in a similar manner to 
the case Am = 1. With an m = 2 perturbation, how- 
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ever, two orbiting fragments develop, contract, and sub- 
sequently encounter a runaway instability in the center 
(bottom right panel in Fig. |8j. Any perturbation with 
different values for Ai and A2 might produce a mixture 
of this spiral arm and binary-system fragmentation in- 
stability. One might argue that a fine-grained parameter 
space study in the space of Aj and A2 is necessary to 
obtain a more complete imderstanding of the remnants. 
However, the reference polytrope is already well inside 
the unstable region of the parameter space. We will see 
in Section fV El that, on a sequence of increasing r/|W|, 
the m = 1 mode dominates first. As has been shown in 
(1, this process of fragmentation is black-hole forming . 
Unfortunately, around the time of horizon formation the 
evolution fails. This appears to be due to the growth of 
small scale features that cannot be adequately resolved 
leading to imphysical oscillations. It seems possible that 
a more robust gauge condition, artificial dissipation (9^ 
or excision techniques [96, 97] could avoid these prob- 
lems. 

To investigate the growth time t of the modes, which, 
in the context of linear theory, is defined by the relation 
Am(f) = A,„{t = 0) exp{t/Tm), we plot the function 
Tm{t) ■- {d\nA/dty^ in the top panel of Fig. [TOl for 
m = 1 and m = 2, using different initial perturbations 
of the reference polytrope as before. As expected for a 
dynamical instability, the growth times are of the order 
of the dynamical timescale. Finally, in the bottom panel 
of Fig.[TO]the mode frequencies cOmit) = d(pm/dt (as ex- 
tracted from the Fourier decomposition of the density 
on a the coordinate radius of highest initial density) are 
plotted versus time. These, and the connected pattern 
speeds w / m, will be important in the discussion of the 
corotation band in Section lV A 31 



2. Some results on the global nature of the instability 

While we do not focus on the global nature of 
the quasi-normal modes of general relativistic quasi- 
toroidal polytropes here, this section will give some in- 
dication about the structure of the instability. Consider 
Fig. [TT] which displays the equatorial distribution of the 
function log^Q{\p{x,y,z, t) — pQ{x,y,z) \ + e), where po is 
the density function of the equilibrium polytrope and 
e = 10^^*^. These logarithmic difference plots exhibit 
the node lines of the unstable mode corresponding to 
m = 1 in the equatorial plane, and show the spiral-arm 
structure of the fragmentation instability. 

The mode amplitudes at different extraction radii in 
the equatorial plane are shown, for a perturbation with 



^ The case = Smi exhibits a "collapse of the lapse" at late times, 
which suggests that in this case a black hole has formed, too. Un- 
fortunately, it was not possible to locate an apparent horizon in this 
case due to numerical difficulties. 




Figure 11: Contour plot of the function logjQ(|jo(x,i/,2, t) — 
PQ{x,y,z) \ + e) in the equatorial plane, where p denotes the 
rest-mass density of the m = 1 evolution shown in Fig. [T] pq 
denotes the the density function of the equilibrium polytrope, 
and e is a small number (e = 10^^" in ttiis plot). 



17 




io"| ^ 1 1 r 




Figure 12: Evolution of the mode amplitudes and Ai in 
the reference polytrope perturbed by A,„ = J„,i, for different 
mode extraction radii. The radius of highest initial density is 
indicated. 

Am = dm\f in Fig. [12] and [13l The development of the 
unstable modes is not very sensitive to the extraction ra- 
dius, at least as long as the amplitude of the dominant 
mode is not close to unity. Fig.[l4]suggests that the local 
mode frequency does not depend strongly on the radius, 
at least within the considerable uncertainties of the plot. 

3. T/ie location of the unstable modes in the corotation hand 

To determine the location of the instability with re- 
spect to the corotation band, we define a coordinate an- 
gular velocity of the initial model by 

Ci{cD) = ocvf - ^f. (18) 

This can be compared to the mode pattern speed 
1/m dcp/dt, which we will assume to be valid for the 
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Figure 13: Same as Fig. [121 but for the mode amplitudes 
and A^. 



whole star (cf . Fig. [14)| , to determine whether a certain 
mode has a corotation point. In Fig. [151 the angular ve- 
locity is plotted in addition to the numerical approxi- 
mation of the location of the m = 1 and m = 2 pattern 
speeds. We find that both modes have corotation points: 
the m = 1 mode near the radius of highest density at 
0.25 Re, and the m = 2 mode near 0.5 .. . 0.6 Re- 



4. Grid resolution and convergence 

For any parameter study with numerical methods, it 
is important to have an understanding of the amount of 
grid resolution needed to extract the physical features 
under consideration. A typical way to gauge this is to 
evolve a system with different resolutions and to com- 
pare the results. For the black-hole forming fragmen- 
tation instability shown in Section fV A li it is expected 
that different phases of the evolution have substantially 
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Figure 14: Local mode frequency tiiy„dij)/dt for the modes 
ni = 1 (top) and m = 2 (bottom) in the reference polytrope, 
extracted at different radii in the equatorial plane. 



different resolution requirements. During the nearly ex- 
ponential growth of the instability at low amplitudes 
(which we will call linear regime), the equilibrium struc- 
ture of the star as a whole needs to be covered appropri- 
ately. The instability is, at first, a low-frequency effect 
on the star, and as such it is not expected to dominate 
the resolution requirements. However, if the fragment 
evolves into a black hole (in the non-linear regime), it 
needs to be resolved with significantly more grid points. 

As explained in Section the star is covered by a 
grid with fixed mesh refinement. Typically, five grid 
patches are used, centered on each other, and with an 
increase of resolution by a factor 2 each. Only the cen- 
tral patch with highest resolution, which covers the re- 
gion of highest density, is 0.75 times as extended as 
the second finest one to reduce artefacts from inter- 
patch boimdaries. To test convergence, the reference 
model has been evolved with 49 x 49 x 25 (grid spac- 
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Figure 15: Angular velocity of the reference polytrope over the 
X axis (black line), and approximate location (with error bar) 
of the pattern speed of the m = 1 mode (upper rectangle), and 
the m = 2 mode (lower rectangle). Both modes are inside the 
corotation band. 



ing h w 0.29M), 65 x 65 x 33 {h w 0.22M), 97 x 97 x 49 
{h w 0.15M) and 129 x 129 x 65 (/z w O.llM) zones per 
outer grid patch; the innermost patch covering the high- 
density toroidal region has 65 x 65 x 33, 97 x 97 x 49, 
129 X 129 X 65 or 193 x 193 x 97 grid zones. Also, the 
initial data was calculated with a grid of 401 x 201, 
601 X 301, 1201 X 601 and 2401 x 1201 zones. 

The results from evolving the reference model with an 
m — 1 perturbation at different resolutions are shown in 
Figs. [TJ [161 HZl and [TS] The convergence of the Hamil- 
tonian constraint has already been discussed in Sec- 
tion IIV Al The time development of the mode ampli- 
tudes has a noticeably different growth only at the low- 
est resolution. The evolution of the maximum of the 
rest-mass density (Fig. [TTt exhibits a similar behaviour. 
Finally, the total rest mass of the system is conserved 
from within 1.4% (lowest resolution) to 0.1% (higher res- 
olutions). The drift in the rest-mass can be explained by 
our use of an artificial atmosphere: The rest-mass den- 
sity of the atmosphere is IQ^^pc = 3.38 ■ 10^^^, which 
corresponds to an approximate total mass of Mo^atmo ~ 
3.8 ■ 10^^ in a domain of coordinate volume 1040'^ (not 
taking into account the volume form). This translates 
into a systematic shift in the total rest mass of the sys- 
tem as apparent when comparing to an evolution with 
a lower atmospheric density (bottom panel in Fig. [TSl l, 
and a drift caused by the intrinsic atmospheric dynam- 
ics and the interaction with the outer boundary. Note 
that this is to be considered a (non-sharp) upper limit 
on the systematic errors induced by the atmosphere, be- 
cause one could always extend the total computational 
domain arbitrarily without affecting the core region sig- 
nificantly. 
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Figure 16: Evolution of the mode amplitudes (top) and A2 
(bottom) for different grid resolutions. The grid sizes in the 
legend refer to the four outer grid patches; the innermost patch 
covering the high-density central toroidal region of the star has 
a resolution of 65 x 65 x 33, 97 x 97 x 49, 129 x 129 x 65 or 
193 X 193 X 97, correspondingly. 



Judging from the resolution study in Figs. [TJ [161 HZI 
and [18l we think that the lowest resolution considered 
here (which already covers the equatorial radius of the 
star with about 60 grid points when comparing to uni- 
form grids) is sufficient to get qualitatively correct re- 
sults. Quantitatively, the errors in rest mass are in a 
range of a few percent. The next highest resolution 
of 65 X 65 X 33 (h w 0.22M) seems accurate to within 
about one percent. This resolution will therefore be 
used for the parameter study below. This is reason- 
able since the structure of the quasi-toroidal models has 
similar features, and therefore similar requirements con- 
cerning resolution. Nevertheless, selected models have 
been tested for convergence independently from the ref- 
erence polytrope. 
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Figure 17: Same as Fig. [161 but for the evolution of the maxi- 
mum of the rest -mass density pmax- 
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Figure 18: Same as Fig.[T6l but for the evolution of the total rest 
mass Mq. The values are normalized to the initial rest mass 
Mq^^ obtained from the initial data solver. The last graph is 
from a simulation with a less dense artificial atmosphere. 



5. Influence of the artificial atmosphere 

The standard artificial atmosphere we employ in our 
simulations has a density several orders of magnitude 
lower than the average density in the star, so we expect 
that it does not influence the dynamical properties of the 
star significantly. The atmospheric density is set in terms 
of the central density of the star: we have used a ratio of 
10^^ in most simulations. To test the influence of this 
parameter, we have evolved the reference model also 
with a ten times lower atmospheric density (i.e. 10~^pc), 
and with an m = 1 perturbation. The results are shown 
in Fig. [TBI and [19l The latter shows that the dominant 
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Figure 19: Evolution of the mode amplitudes Ai and A2 in 
the reference polytrope, for artificial atmospheres of different 
density |0af„,/|0c. 



m = 1 mode is not influenced by the atmospheric set- 
ting, while the m = 2 amplitude shows dependence on 
the atmosphere setting only as long as its amplitude is 
on the level of the numerical noise. 



B. Evolution of the sequence of axis ratios 

The R sequence has been described in Section IIV B 21 
From Table HIl it is apparent that higher values of Vp/Ve 
are connected to lower r/| W|. In Maclaurin spheroids, 
this is related to a stabilization of the initial model. Con- 
sider Fig. |20l The growth time of the modes m = 1 and 
m = 2 increases with lower Vp / r^, which indeed is a sign 
for approaching a limit of stability. 



C. Evolution of the sequence of stiffnesses 

The change in the instability along the G sequence de- 
scribed in Section IIV B 31 is shown in Fig. |21] With in- 
creasing F, and r/|W| decreasing from 0.227 to 0.159 
(cf. Table Unil, both the m = 1 and the m = 2 modes are 
stabilized. The member G2.0 with F = 2 is of special in- 
terest, since this choice is often used to obtain a simple 
polytropic equilibrium model of neutron stars. While 
it is known that strong differential rotation can induce 
bar-mode instabilities in neutron stars II6II |62L |6^ , the 
particular model G2.0 does not appear to be m = 2 un- 
stable (note, however, the limitations of our method to 
determine stability expressed in Section lVGI I. 




Figure 20: Evolution of the mode amplitudes Ai and A2 for 
different members of the R sequence (cf. Table Hit. 



D. Evolution of the sequence of compactnesses 

The mode amplitudes Ai and A2 for different mem- 
bers of the C sequence (cf . Section IIV B 41 and Table IIV|| 
are shown in Fig. |22l The plot demonstrates that differ- 
ent choices of pc do not have a significant effect on the 
growth time of the mode, which is in contrast to the ef- 
fects of F and fp / discussed above. However, while 
the linear development of the mode is similar for dif- 
ferent compactnesses, the non-linear behaviour is not. 
Consider Fig. |23l The reference polytrope CI and the 
model C2, which is about half as compact, both show 
an unbounded growth in the maximum of the density 
and a collapse of the lapse, indicating black hole forma- 
tion. The models C4 and C8, however, appear to avoid 
black hole formation and re-expand after a state of max- 
imum compression. Fig. |24] shows these different types 
of evolution for the model C8 in more detail: In con- 
trast to the black hole-forming case CI, the fragment re- 
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Figure 21: Evolution of the mode amplitudes A\ and Ai for 
different members of the G sequence (cf. Table Ullt. 

expands after the collapse in this evolution. This is not 
unexpected, since the Newtonian limit, for an equation 
of state r = 4 /3, admits a stable equilibrium state of the 
fragment if it has sufficient rotation. We conclude there- 
fore that, even when the growth time of the instability is 
quite similar for stars of different compactness, the out- 
come of the fragmentation can differ drastically.^ 

E. Evolution of quasi-toroidal models of constant central 
rest-mass density 

The structure of the parameter space plane F = 4/3 
and jOc = 10"'' has been discussed in Section flV CI As 
already noted, the necessity of investigating only one 



These results have been confirmed with lower and higher resolu- 
tions. 



Figure 22: Evolution of the mode amplitudes A\ and A2 for 
different members of the C sequence (cf. TablellVt. 



plane is determined primarily by the computational cost 
of three-dimensional general relativistic hydrodynami- 
cal simulations. Also, the choice of the central density 
does not seem to affect the almost exponential develop- 
ment of a non-axisymmetric unstable mode in the lin- 
ear regime considerably, even for very compact quasi- 
toroidal polytropes (cf. Section IV D|l , which is in con- 
trast to axisymmetric modes. Finally, we note that the 
models with pc = 10^^ are already quite compact, with 
Re/M w 10 . . . 100, and rp/M w 2 .. . 70. 

To investigate the stability of these models, the ini- 
tial data indicated by circles in the bottom right panel of 
Fig.|4]have been evolved, imposing a perturbation of the 
form given by eqn. [T7]with A,„ = 1, and with a resolu- 
tion of 65 X 65 X 33 grid points in the outer patches, and 
97 X 97 X 49 in the innermost patch. Selected models 
have been tested against individual m — j perturbations 
with Am = with different resolutions, and different 
densities of the artificial atmosphere, to test consistency 
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Figure 23: Evolution of the maximum of the rest-mass density 
p and the minimum of the lapse function a for different mem- 
bers of the C sequence. 

and convergence. Also, central rest-mass densities dif- 
ferent from lO^"^ were investigated in a few models. 

Fig. |25] gives an overview of the stability properties of 
the selected models. The Latin numbers "I" to "III" refer 
to the highest m with an unstable mode, i.e. in addition 
to the reference polytrope, which belongs to the class 
"II," we find models which are unstable to an m = 3 per- 
turbation, and models which appear to be stable against 
m = 2 (within the restrictions illustrated in Section lVGI I. 
The models denoted with an "A" have been found to be 
unstable to an axisjrmmetric mode, and collapse before 
any non-axisymmetric instability develops. Finally, the 
models marked with "(I)" are either stable or long-term 
unstable with a growth time t ^ t^yn- Each model has 
been evolved for up to 10 f^fy^ to determine its stability. 
This limit is arbitrary, but imposed by the significant re- 
source requirements of these simulations. If no mode 
amplitude exceeds the level of the m = 4 noise during 




Figure 24: Evolution of density in the equatorial plane of the 
model C8. Shown is the decadic logarithm of the rest-mass 
density. The snapshots were taken at times t/t^y„ = (top), 
6.28 (middle) and 8.28 (bottom). In contrast to the more com- 
pact model CI, which is the reference polytrope investigated 
earlier (cf . Fig. |5), the fragment re-expands after a maximal 
compression. 
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Figure 25: Stability of quasi-toroidal models with pc = 10^'' 
(cf. Fig. ID. A Latin number denotes the highest azimuthal or- 
der of the unstable modes, i.e. "1" for m = \ unstable, "11" 
for m = 1,2 unstable, and "111" for m = 1,2,3 unstable. Mod- 
els denoted by "(1)" are either long-term unstable with growth 
times T ^ trfy„, or stable (see text), and models denoted by "A" 
exhibit an axisymmetric instability. The red line in the lower 
left is the approximate location of the sequence // IVfi = 1 (cf. 
Fig. 2), and the three blue lines inside the quasi-toroidal region 
are the approximate locations of sequences with T/\W\ = 0.14 
(right), T/\W\= 0.18 (middle) and T/\W\= 0.26 (left). 



this time, the model is marked with a "(I)." This does 
not imply that the model is actually stable, and we will 
investigate a specific model denoted by "(I)" later. We 
will find it to be unstable to an ffi = 1 mode with slow 
growth (Section lVGI I. 

The additional lines in Fig. |25] are approximate iso- 
lines of the fimctions T/|W| for the values 0.14, 0.18 
and 0.26 and of the function } / for the value 1. 
As long as the models do not rotate too differentially, 
T / 1 W| still seems to be a reasonable indicator of the non- 
axisymmetric stability of the polytropes, even though 
they are quasi-toroidal and relativistic. 

The nature of the non-linear behaviour of models ex- 
hibiting a non-axisymmetric mstability is indicated m 
Fig. |26l We use the evolution of the minimum of the 
lapse fimction to classify the models, see also Fig. |23l 
Models denoted by "B" have a global minimum in the 
lapse, while models denoted by "C" do not. Given that 
the compactness of the models increases with smaller 
axis ratios in this plot (cf. Fig. |4]l, we expect that a 
black hole forms for each member of the "C" class. To 
determine this uniquely, each of these models should 
be tested using the adaptive mesh-refinement technique 
presented in [2], which is, however, beyond the scope of 
this study.^ 



Figure 26: Remnants of the models from Fig. |25l which are 
unstable with respect to non-axisymmetric modes. The non- 
linear behaviour has been analyzed by observing the evolution 
of the function mina (see also Fig.|23l. Models which show a 
minimum in this function are marked by "B" for "bounce," 
while models exhibiting an exponential collapse of the lapse 
are marked by "C" for "collapse." 
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Figure 27: Mode amplitudes in the model A0.1R0.15 (cf. Ta- 
ble E), extracted at the radius of highest initial rest-mass den- 
sity. 



In Fig.|27|to|3TJ we have plotted the mode amplitudes 
Am for selected models (cf. Table |V|. The evolution of 
the model A0.3R0.15 (Fig. |29]l is quite similar to that of 
the reference polytrope. Model A0.6R0.15 is further in- 
side the unstable region, and exhibits also an 7« = 3 in- 
stability: the density evolution of this mode is plotted in 
Fig.|32l The models A0.1R0.50 and A0.3R0.50 are stable 



' We also note that models denoted by "B" might actually form a 



black hole by delayed coUapse. 
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Figure 28: Mode amplitudes in the model A0.1R0.50 (cf. Ta- 
ble |V), extracted at the radius of highest initial rest-mass den- 
sity. 



Figure 30: Mode amplitudes in the model A0.3R0.50 (cf. Ta- 
ble |V), extracted at the radius of highest initial rest-mass den- 
sity. 
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Figure 29: Mode amplitudes in the model A0.3R0.15 (cf. Ta- 
ble |V), extracted at the radius of highest initial rest-mass den- 
sity. 



Figure 31: Mode amplitudes in the model A0.6R0.15 (cf. Ta- 
ble |V), extracted at the radius of highest initial rest-mass den- 
sity. 



within the numerical restrictions mentioned above. 

The model A0.1R0.15 seems to have an unusual evo- 
lution of the mode amplitudes; this also applies to the 
model A0.1R0.20, which is not shown here. The den- 
sity evolution of A0.1R0.15 (Fig. |33ll shows that the 
model has encountered an axisymmetric instability be- 
fore any non-axisymmetric modes can grow to appre- 
ciable amplitudes. We note that both models A0.1R0.15 
and A0.1R0.20 have }/Iyf < 1, in contrast to most 
other models in the parameter space plane considered 
here; only the model A0.1R0.25 has = 0.961. In 

Fig- ESI the isoline }/M^ = 1 is marked. It approxi- 



mately separates the region of axisymmetric from that 
of non-axisymmetric instability. 



F. The location of the instability in the corotation band 

Fig. |34] shows the location of the imstable modes in 
the corotation band, for three different models on the se- 
quence F = 4/3, pc = 10^'^, and A = 0.3. All modes are 
in corotation, but there is evidence that, for decreasing 
T/|W|, the corotation point moves towards the axis of 
rotation. This gives support to the arguments presented 
in [3], where the existence of low-T/| W| and spiral-arm 
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Figure 33: Equatorial density evolution of the model 
A0.1R0.15. Shown is the decadic logarithm of the rest-mass 
density. The snapshots were takes at times t/t^y„ = (top) 
and 6.28 (bottom). The model exhibits an axisymmetric insta- 
bility 



Figure 32: Equatorial density evolution of the model 
A0.6R0.15, with an m = 3 perturbation. (Note that, in Fig.|3T1 
a perturbation with Am = 1 has been used instead). Shown 
is the decadic logarithm of the rest-mass density. The snap- 
shots were takes at times t/ tjy„ = (top), 6.28 (middle) and 
7.60 (bottom). Three fragments develop and subsequently en- 
coimter collapse similar to the two-fragment case (cf. Fig. |8). 
The evolution of the model perturbed with m = 1 and m = 2 
is similar to the corresponding one in the reference polytrope. 



instabilities in differentially rotating polytropes are con- 
nected to corotation^'^. Limitations of resources did not 
permit us to investigate the boundary of the corotation 
region, where growth times of many f ^^^,, ar e expected. 
However, with the results of Sections IV Dl and IV Hi a 
purely Newtonian investigation could be sufficient to 
reproduce the linear regime of the instability. 



There is evidence that, on a sequence of increasing rotation param- 
eter, some modes in the discrete spectrum become unstable when 
entering the corotation band (which has a continuous spectrum), or 
might merge with other modes inside the band and become unsta- 
ble. These are mechanisms not present in uniformly rotating poly- 
tropes. 
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Figure 34: Angular velocity of the polytropes A0.3R0.15 (top), 
A0.3R0.30 (niiddle) and A0.3R0.40 (bottom) over the x axis 
(black line), and approximate location of the pattern speed of 
the m = 1 mode (upper rectangle) and the m = 2 mode (lower 
rectangle), cf. also Fig. [151 
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Figure 35: Long-term evolution of the mode amplitudes for the 
model A0.2K0.45, which is unstable to an m = 1 perturbation. 
The mode, however, grows rather slowly over a time of 20 t^^,,. 



G. Evolution of a model with a slow growth of the m = 1 
instability 

As already discussed, the nature of the class "(I)" 
models in Fig.|25]could not be investigated in detail due 
to the high computational cost when evolving general 
relativistic, three-dimensional models. However, to il- 
lustrate the behaviour in one specific case, a long-term 
simulation of the model A0.2R0.45 has been performed 
(Fig.l35]|. A slowly growing m = 1 instability is appar- 
ent in the evolution, which saturates at high amplitudes 
only after 20 f^y,,. While the m = 1 mode is clearly dom- 
inant, the m = 2 might be unstable as well. A detailed 
investigation of these sequences should be attempted in 
the limit of vanishing compactness, with a Newtonian 
model and preferably, with a cylindrical grid (see also 
the discussion in Section rvTl l. 



H. Evolution of a sequence of models with different 
compactness starting from the boundary between the 
regions "I" and "(I)" 



In Section rVDi we have already studied the influence 
of the compactness on the development of the instabil- 
ity in the reference polytrope. According to the results of 
Section rVEi the reference model is located inside region 
"11" of the parameter space plane for pc = 3.38 ■ 10"^. 
Thus, it is instructive to investigate the effect of com- 
pactness on a model's evolution which is located close 
to the boundary between regions "I" and "(I)" in Fig.|25l 
(although this boundary is not sharply defined). A se- 
lected model, A0.2R0.40, has been extended to the L se- 
quence of constant T, A, and r/|W| (cf. Table IVH i. The 
influence of compactness on the m = 1 mode is illus- 
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Figure 36: Amplitude of the m = 1 mode for different models 
from a sequence with constant T/\\N\ limiting in the model 
A0.2R0.40, which has pc = 10"^. Cf. also Table|Vl] 

trated in Fig. |36l The most compact models Ll and L2 
show a growth of the non-axisymmetric mode already 
early on, but collapse due to an axisymmetric instability 
(both models have } / < 1). The growth rate of the 
non-axisymmetric instability is not very sensitive to the 
compactness, which confirms our findings for models of 
the C sequence. One might therefore be reasonably op- 
timistic that the non-axisymmetric stability properties of 
quasi-toroidal N = 3 polytropes are well-represented 
by Fig. |25l even for a different choice of central rest- 
mass density. The axisymmetric stability, and the ques- 
tion whether the collapse of the fragment will be halted 
or not, is sensitive to pc- 

VI. DISCUSSION 

In this paper we have presented an extension of our 
earlier work on the fragmentation of general relativis- 
tic quasi- toroidal polytropes and the production of black 
holes in this scenario [2]. The central focus was to gain 
an understanding of how various parameters determin- 
ing the structure of the equilibrium polytrope affect the 
development of the instability observed in [2J, and the 
nature of its remnant. In addition, we have investigated 
the location of the unstable modes in the corotation band 
of the differentially rotating models. 

All investigations have been performed using three- 
dimensional numerical simulations in general relativity, 
and assuming the stars to be self-gravitating perfect flu- 
ids with an adiabatic coefficient equal to the polytropic 
constant F. The equations of general relativistic hy- 
drodynamics have been evolved using high-resolution 
shock-capturing methods, and the NOK-BSSN formal- 
ism has been used for the metric evolution. All grids 



use fixed mesh refinement, and impose an equatorial 
plane symmetry. The development of unstable modes 
has been followed by the use of a discrete Fourier trans- 
form of the rest-mass density computed at certain coor- 
dinate radii in the equatorial plane, with a preference on 
the radius of initial highest density. 

The central results are represented in Fig. |25l |26l and 
|34l For a plane of constant rest-mass density pc = 10^'' 
and F = 4/3, we have determined the region where 
quasi-toroidal models become dynamically unstable to 
non-axisymmetric fragmentation. From the structure of 
the space of initial models presented in Fig.Sl it appears 
that there is a rough relation between r/|W| and the 
highest order of unstable modes, at least as long as the 
degree of differential rotation is not too high. Since the 
numerical method is not well-suited to follow the devel- 
opment of instabilities with growth times much longer 
than a dynamical timescale, we could not determine the 
fate of models from class "(I)" with certainty. However, 
we have shown in one specific case that the model is 
actually unstable. In the same manner, a model from 
class "I" could be subject to a slowly growing mode with 
m > 1; however, in this case the m = 1 mode would 
clearly be dominant. 

From the investigation of a sequence emanating from 
the model used in the publication f?], we have found 
that the central rest-mass density pc, which controls the 
compactness of the polytrope, does not affect the al- 
most exponential development of the non-axisymmetric 
instability significantly. This is related to the fact that 
T/|W| is insensitive to pc- However, pc determines the 
nature of the final remnant: While the model in 1 2] forms 
a black hole, two models having one fourth and one 
eighth as much compactness show a re-expansion of the 
fragment after maximal contraction. 

The regions of models in the plane pc = and 
F = 4/3 where such a re-expansion was observed is 
indicated in Fig. |26] by a "B." If one assumes that the 
models not exhibiting this behaviour, marked with "C" 
in Fig. |26l are forming black holes in the same manner 
as shown in [2], then we can draw the following tenta- 
tive picture of black hole formation by fragmentation of 
single stars. 

The nature of the final system, either an almost un- 
perturbed axisymmetric star, a single central black hole, 
single or multiple non-central black holes with a disk, 
or one or several expanding remnants without trapped 
surfaces, depends on the symmetry properties of the 
perturbation, and the location of the equilibrium star 
with respect to three types of surfaces in the space of 
parameters: 

1. Surfaces indicating the onset of the instability of a 
mode of a certain order m. These surfaces might 
be close to isosurfaces of r/|W| as indicated in 
Fig. |25l but the resource requirements of perform- 
ing three-dimensional general relativistic simula- 
tions limit our ability to identify slowly growing 
modes. However, if only modes growing on a 
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dynamical timescale are considered, then r/|W| 
might yield a reasonable indicator of the location 
of the limit surfaces. The compactness of the ini- 
tial model seems to have no significant effect on 
the location of these surfaces, at least for T = 4/3. 

2. Surfaces indicating the onset of an axisymnietric insta- 
bility. In our samples, all models with }/M^ < 
1 were unstable to quickly growing axisymmet- 
ric modes, and hence will likely evolve to central 
black holes. The surface J/M^ = 1 can therefore 
be used as an approximate separator between ax- 
isynmietric collapse and stability (cf. (98*] for a 
more detailed discussion of this point). 

3. Surfaces separating prompt black-hole formation from 
re-expansion. In Fig. |26l an approximate determi- 
nation of such a surface has been attempted. In 
a first approach, and assuming that results for the 
stability of slo wly and i miformly rotating relativis- 
tic polytropes 1991 11001 11011 can be applied to the 
fragments, we expect a fragment with a higher 
compactness and lower rotation rate to be destabi- 
lized, and, given that the geometric development 
of the fragmentation process is similar for different 
choices of compactness of the equilibrium poly- 
trope, that there is a close connection to isosurfaces 
of Re/M and J/M^. However, when comparing 
the structure oi Re/M and //M^ (cf. Fig. HJ and 
|26l the situation appears more complicated, and 
deserves further attention. 

With respect to the question of whether multiple black 
holes may form or not, two comments are in order: First, 
in the unstable systems of class "11" and "III," all growth 
times of modes with different m are of comparable mag- 
nitude, i.e. the nature of the actual time development 
in a specific star will depend on the nature of the per- 
turbation, as already mentioned. On a sequence of in- 
creasing r/|W|, we have found the m = 1 perturba- 
tion to be dominant before higher order modes become 
unstable, which suggest the (off-center) formation of a 
single black hole with a massive accretion disk. How- 
ever, as discussed in the main text, this conclusion ap- 
plies to instabilities with growth times in the order of a 
dynamical timescale. More specifically, our simulations 
do not exclude, on a sequence of increasing r/|W|, a 
slowly growing higher order instability to be effective 
before the m = 1 mode is dynamically unstable. Sec- 
ond, when two fragments are forming and collapsing, 
in all cases considered here a runaway instability devel- 
ops and leads to a central collapse (Fig. jSjl. In this case, 
the gravitational wave signal is expected to resemble the 
ring-down phase of a highly deformed black hole. For 
the m = 1 collapse, we have published an approximate 
prediction for a partial waveform in |2], which suggests 
an amplitude and frequency well inside the LISA sensi- 
tivity if the source is a collapsing supermassive star. 



Concerning the nature of the non-axisymmetric 
mode, we have collected evidence that, along a se- 
quence with decreasing T/|W|, the corotation point 
moves towards the axis of rotation of the polytrope. This 
gives support to the arguments by Watts et al. [3]. To 
also investigate the cases with large growth times, how- 
ever, a Newtonian model, preferably on a cylindrical 
grid, would be of advantage to obtain more detailed re- 
sults i^. 

A final comment is in order concerning the black 
hole remnant. Since the normalized angular momen- 
tum / / of the initial model is greater than unity, the 
resulting black hole, unless it is ejected from its shell, 
may very well be rapidly rotating, spun up by accretion 
of the material remaining outside the initial location of 
the trapped surface. 

Possible future work on this problem can be roughly 
divided into four approaches. First, the nature of the 
low-r/|W| and m — 1 instabilities in quasi- toroidal 
polytropes could be investigated in Newtonian gravity, 
or perhaps using some perturbative approximation of 
general relativity, to determine the location of the corre- 
sponding surfaces in parameter space, and to suggest re- 
gions where the quantity r/| W| is still a good indicator 
for the degree of instability. Since the Newtonian poly- 
tropes can be considered to be limit points of relativistic 
sequences with vanishing compactness, the systematic 
effects of general relativity on their stability properties 
can be determined separately. 

Second, the location of the surfaces separating black 
hole formation from "bounce" behaviour, and its rela- 
tion to the initial compactness Re/ M and the normal- 
ized angular momentum } / M^, needs to be determined 
with more detail, specifically also for different equa- 
tions of state and rotation laws. Could a newly formed, 
rapidly and differentially rotating neutron star fragment 
in this way? We have found no example of this kind 
here, but such a question deserves further attention. A 
more detailed investigation of the "bouncing" models 
might also be interesting, since they could have a rich 
phenomenology, including a possible delayed collapse 
to black holes and complex gravitational wave signa- 
tures. 

Third, a study of the evolution of the black hole and 
its shell would shed light on a number of interesting as- 
pects: (i) The angular momentum of the remnant black 
hole, which is related to the structure of the accretion 
disk and its ability to power jets, (ii) the kick velocity 
of the black hole, and the question whether it may be 
ejected from its host galaxy in some cases, and (iii) the 
gravitational wave signal from the formation process at 
^ + . While we were not able to evolve the black hole 
for an extended time after its formation, recent advances 
in discrete techniques may be able to solve this prob- 
lem, either by adding a sufficient amount of artificial 
dissipation [95] or by employing multi-block grids with 
summation-by-parts operators I961 I97I1 . 

Fourth, to connect more closely to certain astrophysi- 



29 



cal systems, a detailed model of the micro-physical pro- 
cesses, particle transport, and magnetic fields is nec- 
essary in many cases to obtain specific answers. The 
most important bulk property appears to be a change in 
r with density, since this would modify the non-linear 
evolution of the fragmentation significantly. In the spe- 
cific case of core collapse, results in this context have 
been obtained already li48i.,63il . 
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